!****************************************************************************************
SUBROUTINE IRFS3
! compute responses in initial period

USE GLOBAL

IMPLICIT NONE

INTEGER::nirf
PARAMETER (nirf=2)

REAL(8),DIMENSION(nirf,1)::vt,at,pt,ht,ctt,gnt,bt,acct,deft,rt,wt,taut,bort,rtc,byt
REAL(8),DIMENSION(nirf,1)::yt,ynt,ytt,cnt,ct,gaint

REAL(8)::vt0,at0,gnt0,ctt0,cnt0,c0,util0,util0ss,ht0,pt0,wt0,bt0,taut0,rt0,rt0c,b00,a_old,a_oldss,a0,a0ss 
REAL(8)::at0ss,gnt0ss,ctt0ss,ht0ss,pt0ss,wt0ss,bt0ss,taut0ss,rt0ss,vt0ss  
REAL(8)::cnt0ss,ct0ss,ynt0ss,ytt0ss,yt0ss,yttv0ss,ytv0ss
REAL(8)::vt2(nt2),at2(nt2),gnt2(nt2),ctt2(nt2),ht2(nt2),cnt2(nt2),pt2(nt2),wt2(nt2)
REAL(8)::bt2(nt2),rt2(nt2),taut2(nt2),yt2(nt2),ytv2(nt2),ynt2(nt2),ytt2(nt2),yttv2(nt2),ct2(nt2),dummy1(nirf),bss

INTEGER::i_at(nirf,1),i_bt(nirf,1),i_at0(nirf,1),i_a0,bort0,acct0,acct02,deft0,bort0ss,acct0ss,deft0ss
INTEGER::i,j,k,t,i_a,i_b,i_ass,i_bss,i_aut0(1),i_at0ss(nirf,1),i_b0(1),i_b0ss,i_at_initial

REAL::eps_a_vec(nirf)
REAL(8)::eps_a,aahat,aa1,utilc,utilg,utilcss,utilgss

REAL(8)::qt0,qt0ss,qt0c,rt0ssc,q_fun,probdef_fun
EXTERNAL::q_fun,probdef_fun

1002 FORMAT(1000F7.2)

seed=1234658753
CALL RNSET (seed)

i_a0 = 1
a0ss = mua/(1d0-rhoa) !start from unconditional mean

!Upload solution with and without spending cuts
CALL INITIAL_COMMIT20

!To compute counterfactual prices, upload baseline value functions and prices
CALL INITIAL_mc

i_bss = i_b0ss      !steady-state debt
bss  = bgrid2(index_global)
b00  = bss 
i_b0 = LOCATE2(bgrid2,b00)

PRINT*, 'Initial (a,b) =',a0ss,b00,bgrid2(index_global),index_global

!i_at_initial = 251 
!For COMMIT_yT
i_at_initial = index_aglobal
i_at0(1,:) = i_at_initial
 
PRINT*,'indices',index_aglobal,index_global
a0ss=agrid2(i_at0(1,1))
i_at0ss = i_at0  
a0 = a0ss 

DO j=1,1

i_b0 = index_global 
bss  = bgrid2(i_b0(1)) 
b00  = bss !initial debt

a_oldss  = a0ss
a_old    = a0 

i_b   = i_b0(1)
i_bss = i_b0(1) 

acct0 = 1   !access in initial period
deft0 = 0
bort0 = 1
    
acct0ss= 1   !access in initial period
deft0ss= 0
bort0ss= 1

CALL r8vec_uniform_01(nirf,seed,dummy1)    !reentry  
CALL RNNOR(nirf,eps_a_vec)

dummy1(1)= 1d0        !no reentry after default

IF ((def1_20(i_b,i_at0(1,j)).EQ.0).OR.(def1(i_at0ss(1,j),i_bss).EQ.0)) THEN
    PRINT*,'Default occurs',j,i_at0(1,j),i_b
    PRINT*,'Default occurs2',i_b,i_at0(1,j),i_at0ss(1,j),i_bss
    PRINT*,'def',def1_20(i_b,i_at0(1,j)),def1(i_at0ss(1,j),i_bss)
    PRINT*,'vaut,vcon_20',vaut1_20(i_at0(1,j)),vcon1_20(i_b, i_at0(1,j))
    PRINT*,'vaut,vcon1',vaut1(i_at0(1,j)),vcon1(i_b,i_at0(1,j))
    GOTO 1234
    PAUSE
    STOP
ENDIF

DO i=1,nirf

    eps_a = 0d0*sigmaa*DBLE(eps_a_vec(i))

    i_a = i_at0(i,j)
    i_ass=i_at0ss(i,j)

    acct02 = acct0
    deft0  = 0
    bort0  = 1
    deft0ss= 0
    bort0ss= 1

    !with spending cut
    IF ((acct0.EQ.0).AND.(dummy1(i)>phi)) THEN
        bort0 = 0 !bort(i)=0
        acct0 = 0 !acct(i+1)=0
    ELSEIF ((acct0.EQ.0).AND.(dummy1(i).LE.phi)) THEN
        bort0 = 1 !bort(i)=0
        acct0 = 1 !acct(i+1)=0 
    ENDIF
    
    IF (acct0.EQ.1) THEN 
        IF (def1_20(i_b,i_a).EQ.1d0) THEN
            bort0 = 1 !bort(i)=1
            acct0 = 1 !acct(i+1)=1
        ELSE
            bort0 = 0  !bort(i)=0
            acct0 = 0  !acct(i+1)=0
            deft0 = 1  !deft(i)=1
        ENDIF
    ENDIF
        
    IF (bort0.EQ.1) THEN    !bort(i)
        vt0  = vcon1_20(i_b, i_a)
        at0  = DEXP(agrid2(i_a))
        gnt0 = gn1_20(i_b,i_a)
        ctt0 = ct1_20(i_b, i_a) 
        ht0  = h1_20(i_b, i_a)
        pt0  = (1d0-omegac)/omegac*(ctt0/(ht0**alpha-gnt0))**(1d0+mu)
        wt0  = alpha*pt0*ht0**(alpha-1d0)
        taut0= taxrate*(at0+in*pt0*ht0**alpha)
        qt0  = q1_20(i_b, i_a)
        rt0  = (1d0+q1_20(i_b, i_a)**(-1d0)-lambda)/(1d0+r)-1d0
        bt0  = b1_20(i_b,i_a)
        i_b  = bp1_20(i_b, i_a)
        !counterfactual prices and spreads
        qt0c = q_fun(bt0,DLOG(at0))
        rt0c = (1d0+qt0c**(-1d0)-lambda)/(1d0+r)-1d0

        acct0 = 1
    ELSEIF (bort0.EQ.0) THEN
        vt0  = vaut1_20(i_a)
        at0  = DEXP(agrid2(i_a))
        gnt0 = gnaut1_20(i_a)
        ctt0 = ctaut1_20(i_a)
        ht0  = haut1_20(i_a)
        pt0  = (1d0-omegac)/omegac*(ctt0/(ht0**alpha-gnt0))**(1d0+mu)
        wt0  = alpha*pt0*ht0**(alpha-1d0)
        taut0= taxrate*(at0+in*pt0*ht0**alpha)
        rt0  = -1d0
        rt0c = -1d0 
        bt0  = 0d0
        i_b  = i_bzero2
    ENDIF
    
    !without spending cut
    IF ((acct0ss.EQ.0).AND.(dummy1(i)>phi)) THEN
        bort0ss = 0 !bort(i)=0
        acct0ss = 0 !acct(i+1)=0
    ELSEIF ((acct0ss.EQ.0).AND.(dummy1(i).LE.phi)) THEN
        bort0ss = 1 !bort(i)=0
        acct0ss = 1 !acct(i+1)=0
    ENDIF
    
    IF (acct0ss.EQ.1) THEN       
        IF (def1(i_ass,i_bss).EQ.1d0) THEN
            bort0ss = 1 !bort(i)=1
            acct0ss = 1 !acct(i+1)=1
        ELSE
            bort0ss = 0  !bort(i)=0
            acct0ss = 0  !acct(i+1)=0
            deft0ss = 1  !deft(i)=1
        ENDIF
    ENDIF

    IF (bort0ss.EQ.1) THEN    !bort(i)
        vt0ss  = vcon1(i_ass,i_bss)
        at0ss  = DEXP(agrid2(i_ass))
        gnt0ss = gn1(i_ass,i_bss)
        ctt0ss = ct1(i_ass,i_bss)
        ht0ss  = h1(i_ass,i_bss)
        pt0ss  = (1d0-omegac)/omegac*(ctt0ss/(ht0ss**alpha-gnt0ss))**(1d0+mu)
        wt0ss  = alpha*pt0ss*ht0ss**(alpha-1d0)
        taut0ss= taxrate*(at0+in*pt0ss*ht0ss**alpha)
     
        bt0ss   = bbp1(i_ass,i_bss)
!        qt0ss   = q_fun(bt0,DLOG(at0))
!        rt0ssc  = (1d0+qt0ss**(-1d0)-lambda)/(1d0+r)-1d0
        !qt0ss   = q_fun(bt0ss,DLOG(at0))
        !rt0ss   = (1d0+qt0ss**(-1d0)-lambda)/(1d0+r)-1d0
        rt0ss   = (1d0+q1(i_ass,i_bss)**(-1d0)-lambda)/(1d0+r)-1d0
        i_bss   = bp1(i_ass,i_bss)
        acct0ss = 1 
    ELSEIF (bort0ss.EQ.0) THEN
        vt0ss  = vaut1(i_ass)
        at0ss  = DEXP(agrid2(i_ass))
        gnt0ss = gnaut1(i_ass)
        ctt0ss = ctaut1(i_ass)
        ht0ss  = haut1(i_ass)
        pt0ss  = (1d0-omegac)/omegac*(ctt0ss/(ht0ss**alpha-gnt0ss))**(1d0+mu)
        wt0ss  = alpha*pt0ss*ht0ss**(alpha-1d0)
        taut0ss= taxrate*(at0+in*pt0ss*ht0ss**alpha)
        rt0ss  = -1d0
        bt0ss  = 0d0
        i_bss  = i_bzero2
    ENDIF
    
    vt(i,j)  = vt0-vt0ss
    at(i,j)  = at0-at0ss
    gnt(i,j) = gnt0-gnt0ss
    ctt(i,j) = ctt0-ctt0ss
    ht(i,j)  = ht0-ht0ss                      
    pt(i,j)  = pt0-pt0ss
    wt(i,j)  = wt0-wt0ss
    i_bt(i,j)= i_b
    bt(i,j)  = bt0 - bt0ss
    byt(i,j) = bt0/(at0 + pt0*ht0**alpha) - bt0ss/(at0ss + pt0ss*ht0ss**alpha) 
    rt(i,j)  = rt0-rt0ss
    rtc(i,j)  = rt0-rt0c
    
    taut(i,j)= taut0-taut0ss
    deft(i,j)= deft0-deft0ss
    acct(i,j)= acct02-acct0ss
    bort(i,j)= DBLE(bort0)-DBLE(bort0ss)
    bt(i,j)  = bt(i,j)/(lambda+r)
    
    ctt(i,j) = 100d0*ctt(i,j)/ctt0ss
    gnt(i,j) = 100d0*gnt(i,j)/gnt0ss
    bt(i,j)  = 100d0*bt(i,j)/(bt0ss/(lambda+r)) 
    pt(i,j)  = 100d0*pt(i,j)/pt0ss    
    
    cnt0ss = ht0ss**alpha-gnt0ss
    IF (mu.EQ.0d0) THEN
        ct0ss = (ctt0ss**(omegac))*(cnt0ss**(1d0-omegac))
    ELSE
        ct0ss = (omegac*ctt0ss**(-mu)+(1d0-omegac)*cnt0ss**(-mu))**(-1d0/mu)
    ENDIF

    IF (alphac.EQ.1d0) THEN
        aahat  = 1d0
        aa1    = 1d0
    ELSE
        aahat  = (ht0ss+(1d0-ht0ss)*alphac)**(-1d0)
        aa1    = ht0ss+(1d0-ht0ss)*(alphac**(1d0-sigma))
    ENDIF

    IF (sigma.EQ.1d0) THEN
        utilcss = (aahat**(1d0-sigma))*aa1*DLOG(MAX(ct0ss,1d-8))
    ELSE    
        utilcss = (aahat**(1d0-sigma))*aa1*MAX(ct0ss,1d-8)**(1d0-sigma)/(1d0-sigma)
    ENDIF

    IF (sigmag.EQ.1d0) THEN
        utilgss = DLOG(MAX(gnt0ss,1d-8))
    ELSE
        utilgss = MAX(gnt0ss,1d-8)**(1d0-sigmag)/(1d0-sigmag)
    ENDIF

    util0ss = (1d0-chi)*utilcss+chi*utilgss 
    
    ynt0ss = ht0ss**alpha
    ytt0ss = at0ss
    yt0ss  = ytt0ss + pt0ss*ynt0ss 
    yttv0ss= ytt0ss         !value-added
    ytv0ss = yt0ss 
    
    cnt0 = ht0**alpha-gnt0
    cnt(i,j) = cnt0-cnt0ss
    cnt(i,j) = 100d0*cnt(i,j)/cnt0ss     

    IF (mu .NE. 0d0) THEN
        c0 = (omegac*ctt0**(-mu)+ (1d0-omegac)*cnt0**(-mu))**(-1d0/mu)
    ELSE
        c0 = (ctt0**(omegac))*(cnt0**(1d0-omegac))
    ENDIF

    IF (alphac.EQ.1d0) THEN
        aahat  = 1d0
        aa1    = 1d0
    ELSE
        aahat  = (ht0+(1d0-ht0)*alphac)**(-1d0)
        aa1    = ht0+(1d0-ht0)*(alphac**(1d0-sigma))
    ENDIF

    
    IF (sigma.EQ.1d0) THEN
        utilc = (aahat**(1d0-sigma))*aa1*DLOG(MAX(c0,1d-8))
    ELSE    
        utilc = (aahat**(1d0-sigma))*aa1*MAX(c0,1d-8)**(1d0-sigma)/(1d0-sigma)
    ENDIF

    IF (sigmag.EQ.1d0) THEN
        utilg = DLOG(MAX(gnt0,1d-8))
    ELSE
        utilg = MAX(gnt0,1d-8)**(1d0-sigmag)/(1d0-sigmag)
    ENDIF

    util0 = (1d0-chi)*utilc+chi*utilg 

    gaint(i,j)= (vt0ss-vt0)/(util0ss-(vt0ss-vt0))

    ynt(i,j) = ht0**alpha -ynt0ss
    ytt(i,j) = at0 -ytt0ss
    yt(i,j)  = ytt(i,j) + pt0*ynt(i,j)-yt0ss 
    
PRINT*,' '
PRINT*,'*************************'
PRINT*,'pn',pt(i,j)
PRINT*,'debt',bt(i,j)
PRINT*,'ct',ctt(i,j)
PRINT*,'unemp',-ht(i,j)
PRINT*,'gN',gnt(i,j)
PRINT*,'spreads',rt(i,j)
PRINT*,'gain',gaint(i,j)
PRINT*,'*************************'

!Storing the results in txt files:
!!For welfare gains from non-state-contingent spending cuts as function of the cut size
!OPEN(101,FILE='output//welfare//welf_gains_state_cutsize_bluntcuts.txt',ACTION='WRITE',POSITION='APPEND')    
!WRITE(1011,'(F18.4,X,F18.3,X,F18.8)') bgrid2(index_global), sizecut, gaint(i,j)
!********************************************************************************************
!!For welfare gains from state-contingent spending cuts as function of the cut size
!OPEN(101,FILE='output//welfare//welf_gains_state_cutsize_statecontingentcuts.txt',ACTION='WRITE',POSITION='APPEND')    
!WRITE(1011,'(F18.4,X,F18.3,X,F18.8)') bgrid2(index_global), sizecut, gaint(i,j)
!********************************************************************************************
!!For welfare gains from non-state-contingent spending cuts as function of current yT
!OPEN(101,FILE='output//cons//welf_gains_blunt_yt.txt',ACTION='WRITE',POSITION='APPEND')    
!WRITE(101,'(F18.4,X,F18.8)') agrid2(index_aglobal), gaint(i,j)
!********************************************************************************************
!!For welfare gains from non-state-contingent spending cuts as function of current debt
!OPEN(101,FILE='output//welfare//welf_gains_blunt_debt.txt',ACTION='WRITE',POSITION='APPEND')    
!WRITE(101,'(F18.4,X,F18.8)') bgrid2(index_global), gaint(i,j)
!********************************************************************************************
!!For welfare gains from state-contingent spending cuts as function of lower bound on applicable total income
!OPEN(101,FILE='output//welfare//welf_gains_lbound.txt',ACTION='WRITE',POSITION='APPEND')    
!WRITE(101,'(F15.3,F15.3,F12.3,F15.8)') bcutmin, bcutmax, gaint(i,j)
!********************************************************************************************
!!For welfare gains from state-contingent spending cuts as function of upper bound on applicable total income
!OPEN(101,FILE='output//welfare//welf_gains_ubound.txt',ACTION='WRITE',POSITION='APPEND')    
!WRITE(101,'(F15.3,F15.3,F12.3,F15.8)') bcutmin, bcutmax, gaint(i,j)
!********************************************************************************************
CLOSE(101)

1234 RETURN


ENDDO  !i
ENDDO  !j
   
PRINT*, ' '
PRINT*, 'IMPUlSE RESPONSES EXPORTED'
PRINT*, ' '
    
END SUBROUTINE IRFS3
